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ABSTRACT 

In this paper we study possible signatures of binary planets or exomoons on 
the Rossiter-McLaughlin (R-M) effect. Our analyses show that the R-M effect 
for a binary planet or exomoon during its complete transit phase can be divided 
into two parts. The first is the conventional one similar to the R-M effect from 
the transit of a single planet, of which the mass and the projected area are 
pLf the combinations of the binary components; and the second is caused by the 

orbital rotation of the binary components, which may add a sine- or linear-mode 
deviation to the stellar radial velocity curve. We find that the latter effect can 
be up to several or several ten ms~^. By doing numerical simulations as well 
as analytical analyses, we illustrate that the distribution and dispersion of the 
latter effects obtained from multiple transit events can be used to constrain the 
dynamical configuration of the binary planet, such as, how the inner orbit of the 
binary planet is inclined to its orbit rotating around the central star. We find 
that the signatures caused by the orbital rotation of the binary components are 
more likely to be revealed if the two components of binary planet have different 
masses and mass densities, especially if the heavy one has a high mass density 
and the light one has a low density. Similar signature on the R-M effect may also 
be revealed in a hierarchical triple star system containing a dark compact binary 
and a tertiary star. 

Subject headings: planetary systems — planets and satellites: detection — plan- 
ets and satellites: dynamical evolution and stability — planets and satellites: 
fundamental parameters — stars: formation — white dwarfs 
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1. Introduction 

Planets are not alone. For example, most of the planets in the Solar System have 
surrounding satellites or moons; not only some asteroids, but also many of the Kuiper Belt 
objects (e.g., Pluto and Charon) are in binaries; and recent observations have also revealed 
miscellaneous exop lanetary worlds, incl u ding multi-planetary systems (e.g., Kepler-11, HP 
10180, Gliese 581; llissauer et aPboili Ihovis et all boiol : IVogt et al.lboill : IPorveille et aL 



20111 ). Although the existence of exomoons or binary exoplanets have not been detected, 
various detection methods have been proposed, such as through transit light curves, transit 
timing and durat ion variations, direct imaging, microlen s ing, or Doppler spectroscopy of the 



host planet (e.g., iKipping et all 120091 : ISimon et al.ll2010l . l2012l : ISartoretti k Schneiderlll999 



and references therein). In this paper, we investigate possible observational signatures that 
a binary exoplanet system (if any) or an exoplanet plus one moon sy stem would have on 
the Rossiter-McLaughlin (R-M) effect ( iMcL au ghlini 1 1 9 241 : lRossiterlll924l ) and how the orbital 
configurations of the system could be revealed through the signatures. Por simplicity, below 
we also call the one exoplanet plus one moon system as a binary planet, although the two 
components have very different masses. Existence of binary planets and statistics on their 
dynamical configurations should shed new light on the formation and evolution of planetary 
systems and the search for a habitable world. 

When a planet transits in front of a rotating star, it blocks part of the light emitted 
from the stellar surface, and the blocked region shifts with the transiting of the planet. As 
different parts of the stellar surface may have different line-of-sight velocities to the observer, 
the shifting of the blocked region results in either blueshift or redshift of observed stellar 
spectral lines and further the deviation of the inferred radial velocity of the stellar motion 
(i.e., the R-M effect) El The deviation of the stellar radial velocity provides a way to measure 
the misalignment between the stellar spin and the planetary orbital angular momentum, 
and recent measurements of the R-M effects have found that some exoplanets a re on highly 



inclined orbits relative to the spiri of the star (e.g., retrograde or polar orbits; I Winn et al. 



2009 



Collier Cameron et al.l 120101 ). Por the similar reasons, if the transiting planet is a 
binary, the binary with different physical properties and orbital configurations (e.g., radii, 
inclination) may block the stellar surface in different ways, and hence the shift of stellar 
spectral lines and the deviation of its radial velocity curve may have different signatures 
from those expected from the transit of one single planet. The gravitation from a binary 
planet is also different from that from one single planet with the same total mass, but the 



^Note that extraction of the velocity devi a-tion from stellar spectrum hne profiles i nvolves some detailed 
techniques in the modeling of the effect (e.g., Hirano et al. 2011 : Albrecht et al. 20071 ) 
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resulted deviation in the dynamical motion of the star is generally small within each transit 
duration (see justification in Section 1331 below) . In this paper, we isolate and illustrate the 
effects on the observed stellar radial velocity due to the different light blocking ways. 

Most of the modeling of the R-M effect for planetary systems were developed for a 
single plan et rotating around a star. The R-M effect for exomoon systems was numerically 



modeled by lSimon et al.l ( l2010l ). where the exomoon has a much smaller mass and radius than 
its host planet and the moon effects are modeled as a small perturbation added to the effect 
caused by a single planet. Their modeling includes the fitting to the numerous dynamical 
parameters of the systems, including instantaneously and fast evolving ones. In this paper 
we present a comprehensive analysis and investigation on the R-M effect for a binary planet 
system, where the satellite mass is not limited to be small but can be comparable to the 
host planet mass. We include the effects from the dynamical evolution of the binary planet 
system during multiple transits and their different dynamical configurations. Our detailed 
treatments average out the effects from some instantaneously changing dynamical angles 
and include the evolutionary patterns of those relatively fast changing ones, so that we can 
focus on the effects from different inner orbital inclinations of the binary planet relative to 
its orbit rotating around the central star. 

The paper is organized as follows. In Section [2l we describe the geometric configuration 
of the system to be studied in this paper (i.e., a binary planet transiting in front of a star) 
and related dynamical approximations. In Section [31 we investigate how the R-M effect 
is affected by a binary planet and how the different orbital configurations could be inferred 
from the deviation of the stellar radial velocity curves, together with the transit of the stellar 
light curves. We explore the parameter space of binary planets that are likely to be revealed 
in observations. We also extend the results to hierarchical triple star systems in which a 
binary star (e.g., a compact stellar remnant plus a brown dwarf or planet) is transiting in 
front of a tertiary star gravitationally bound to the system. Section H] contains a summary 
and a discussion. 



2. Geometric configuration and dynamical description of the system 

Consider that a binary planet is rotating around a star (see Fig. [1]). For convenience, 
we shall call the binary planet as the "inner" binary, and call the star and the center of the 
mass of the binary planet as the "outer" binary. We denote the mass and the radius of the 
star as and i?*. The projected area of the star onto the sky is A^, = nRl, and it has 
a surface brightness of The angular velocity of the stellar spin is 1)*. We denote the 
component masses of the inner binary as mi and m2 (mi > m,2), and the component radii 
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as Ri and R2, respectively. Here we have mi + m2 ^ m*. The projected areas of the two 
components onto the sky are Ai = 'kR\ and A2 = nR"^, respectively (see the parameter list 
in Table HD. 

We describe the dynamical motion of the system through the two components: (i) the 
orbital motion of the outer binary, and (ii) the orbital motion of the inner binary. The orbital 
motion of the inner binary includes the precession of its orbital angular momentum around 
the orbital angular momentum of the outer binary, as described below. 

The orbital configuration of the system is indicated in a reference frame as shown in 
Figure [TJ In Figured! the center of the star is located at the origin O. The y-axis is directed 
toward the observer, and the z-axis is chosen so that the stellar spin axis lies on the y-z plane. 
The inclination angle of the stellar spin relative to the y-axis is denoted by (0 < < tt). 
The unit vector of the orbital angular momentum of the outer binary is denoted by n, and 
we define A by the angle between the z-axis and the projected vector of n onto the x-z plane 
(0 < A < 2tt). The orbital inclination angle to the observer i is defined by the angle between 
n and the y-axis. Thus, we have n = (sin A sin i, cos i, cos A sin i). The outer binary has a 
semimajor axis of a and an angular velocity to = [G{m^ + mi + m2)/a^Y^'^. For simplicity, 
the eccentricity of the outer binary e is assumed to be zero, unless otherwise specified. 

Similarly, we denote the unit vector of the orbital angular momentum of the inner binary 
by n', and we define A' by the angle between the z-axis and the projected vector of n' onto 
the x-z plane, and define the orbital inclination angle i' to the observer by the angle between 
•n' and the y-axis. We have n' = (sin A' sin cos i', cos A' sin z'). The angle between n and n' 
is denoted by 

6 = arccos(n ■ n'), < < vr. (1) 

The semimajor axis and angular velocity of the inner binary are denoted by d and u' = 
[G{mi + m2) / d^Y^'^ , respectively. The eccentricity of the inner binary e' is assumed to be 
zero, unless specially discussed in some cases below. Note that the semimajor axis of the 
inner binary is limited by the Hill radius dn and the Roche limit d^, i.e., 

" ^ " (^) " - im) (^) ( t) " 

/ \ 1/3 

d>dn^ Ri h^] , (3) 

where Mq is the solar mass, Mj is the Jupiter mass, pi (= i = 1,2) is the mass density 
of each component of the inner binary. An inner binary with a larger semimajor axis will be 
tidally broken up by the gravitation from the star. 
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Table 1. List of parameters. 



Object 


Property 


O 11 

Symbol 


Central star 


mass 






radius 






sky-pro iected area 






surface briffhtnpss 






unit vector of spin 






inclination of spin 






spin angular velocity 




Binary planet 


mass 






radius 


llA . ±1,0 




sky-projected area 


^1, ^2 




mass densitv 


On 


Orbit of outer binary 


semimajor axis 


d 




eccentricitv 


g 




unit vector of angular momentum 


— * 

n 




an(?ular velocitv 






sky-projected angle between n and 


A 




inclination to the line of sight 


i 


Orbit of inner binary 


semimajor axis 


d 




distance to center of mass 


di, d2 




eccentricity 


e' 




unit vector of angular momentum 


n' 




angular velocity 


u' 




sky-projected angle between n' and 


X' 




inclination to the line of sight 


i' 


Precession of n' around n 


angle between n and n' 


e 




angular velocity 





Note. — 



See details in Section [2] and Fig. [H 
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Fig. 1. — Schematic diagram for the orbital configuration of the system in which a binary 
planet is rotating around a star. The star is located at the origin O. The two components 
of the binary planet are denoted by "1" and "2" , and their center of mass is labeled by O' . 
For brevity, the star and the binary planet are all indicated by points, without illustrating 
their physical radii in the figure. The observer is along the direction of the |/-axis. The 
spin direction of the star n* is located at the y-z plane. The round curve centering at O 
represents the orbit of the center of mass of the binary planet rotating around the star. The 
n represents the unit vector of the orbital angular momentum of the outer binary, and the 
n' represents the unit vector of the orbital angular momentum of the inner binary (binary 
planet). The angle between n and n' is defined by 6 in Tabled! Point P is the pericenter 
of the orbital plane of the outer binary if the orbit is eccentric, and point Q is located on 
the intersection line of that orbital plane and the plane where both the vector n and the 
7/-axis are located. The angle between and {ZPOQ not labeled in the figure for 
visual clarity) is used in equation ( 1T5]1 later. For other labeled angles, see also Table [T] and 
Section [2j 
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In general, the orbital parameters of the outer binary and the inner binary may evolve 
with time under the three-body interactions of the star and the binary planet. But under 
some conditions, their dynamical motion can be much simplified as follows. 

• The orbital parameters of the outer binary (a, e, n, i, A) can be approximately constant, 
if the semimajor axis of the inner binary is much smaller than the Hill radius dn- 

• The semimajor- axis of the inner binary d can also be approximately constant if it is 
much smaller than the Hill radius. Both the inclination i' and the projected angle A' 
of the inner binary may change with time due to the evolution of n'. The evolution of 
n' depends on the angle 6 betwee n the orbital ang ular momenta of the inner and the 



outer binaries (i.e., n' and n; e.g., iFord et al.ll2000l ). 



(a) If the angle 6 or 180° — 6^ is small (e.g., < 40°, the Kozai angle; lKozailll962l ). the n' 
precesses around n approximately at a constant angular velocity Q ~ u^cosO/u', 
and the angle 6 can also be approximately constant. The angles i' and A' change 
periodically with the precession of n'. As one transit duration St (~ uj~^R^/a) is 
generally much shorter than the precession timescale of n', the n' (and i'. A') can 
be approximated as constant during each transit duration. 

(b) If the angle 9 is about in the range from 40° to 140°, the Kozai mechanism 
affects the evolution of the orbital parameters of the inner binary, in which the 
angle 9 and the eccentricity e' exchange at a period ~ 21100' /u'^ due to angular 
momentum transfer and the conservation of the quantity Ck = Vl — e'^ cos 9 (see 
an example shown in Fig. H] below). During the oscillation, the eccentricity of the 
inner binary can be induced to high values close to 1, and thus the inner binary 
is likely to be destroyed by collision of its two components. If there exist other 
moons in the system, one component of the binary is also likely to be ejected 
from the system or collide with one of the moons. The Jovian system is such 
an example influenced by the Kozai mechanism, and almost all the inclinations 
of their moons are out of the angle range. Here we ignore other effects from the 
planet (e.g., tides, the general relativistic effect), which would limit the influence 
of the Kozai mechanism. The limitation would become relatively significant for a 
system with small d but large a. For example, Uranus is farther from the sun and 
its mner moons (with -R*/a ^ 2.4 x 10""^ and a/d 5x 10^ — 6 x 10^) are on polar 



orbit s relative to its orbital plane surrounding the sun (e.g., iMurray &: Dermott 
1999h . 



Below we do not focus on case (b). 
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3. Transit of a binary planet and its signatures on the R-M effect 



3.1. Transit light curves 



When a binary planet transits in front of a star (see a schematic diagram in Fig. [2]), the 
stellar light curves may be imprinted with signatures of the binary planet. We illustrate one 
example of the normalized transit light curve in Figure |3] (see panel a) . The observational 
transit curve is obtained by L = j J I{x, z)dxdz, where I{x, z) is the observational stellar 
surface brightness at position (x, z). During the transit, part of the stellar surface is blocked 
by the planets, and we set 



For simplicity, the limb-darkening effect of the stellar surface brightness is ignored in this 
paper. We use a full three-body simulation to obtain the dynamical motion of the system. 
The dynamical motion of the binary planet relative to the star determines the shifting of the 
blocked region during the transit. The following phases during the transit are illustrated in 
Figures [2] and [31 

• Ingress phase: at the beginning of the transit, at least one component of the binary 
planet starts to block the stellar light, but the projected areas of the two components 
have not been fully enclosed by the projected stellar surface (e.g., the "AC" part in 



• Complete transit phase of the binary planet: both the projected areas of the two 
components have been fully enclosed by the projected stellar surface (e.g., the "CE" 
part in Fig. [3]). Generally the two planets are more likely to spend most of the transit 
time in that phase when i is close to 7r/2. The d < 2i?^, can be roughly taken as a 
condition for the occurrence of this phase during the transit. If the binary planet has a 
sufficiently large a;', it is likely that during the transit, their orbital evolution leads to 
the evolution of their projected areas from non-overlap to overlap, and to non-overlap 
again, for which a "bulge" (the "DD'" part in Fig. [3]) are shown in the light curve. 

• Egress phase: at least one component of the binary planet has transited to the other 
end of the projected stellar surface and its projected area is not fully enclosed by the 
projected stellar surface again (e.g., the "EF" part in Fig. [3]). The "E'F'" fiat part 
of the light curve in Fig. [3^ represents the period in which one component has fully 
moved out of the projected stellar surface, but the other one is still completely inside. 




unblocked region with 
0, otherwise. 



+ z- 



(4) 



Fig. ED- 
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e f g h 



Fig. 2. — Panels (a)-(g): schematic diagram for transit of a binary planet (small shaded 
circles) in front of a star (big open circle) . The observer is located in the direction of pointing 
out of the paper surface. The binary planet is transiting from left to right. Panels (a)-(c): 
ingress phase; panel (d): complete transit phase of the binary planet; panels (e)-(g): egress 
phase. In this example, the rotation of the binary planet is fast enough so that the small 
component has changed its relative position from the left side to the right side of the big 
one during the complete transit phase (d). See details in Section IXTl Panel (h): schematic 
diagram for the geometry of the system during the complete transit phase of a binary planet. 
The sizes of the star and the binary planet are all illustrated in the panel. The dashed lines 
represent the projection of the labeled vectors onto the sky. See also Table [1] for the meaning 
of labeled parameters. 
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Fig. 3. — Example of transit light curve and stellar radial velocity anomaly for transit of a 
binary planet in front of a star. The reflected light from the surface of the binary planet is 
ignored for simplicity. Panel (a): normalized observational stellar light curve L{t)/{TTl^,Rl). 
Panel (b): observational radial velocity anomaly of the star. The left y-axis in panel (b) 
represents the dimensionless velocity values, and the right y-axis represents the physical 
values of the Af* by assuming Q^R^smi^, = 5km s~^ (similarly for Avf, in Fig. [S] and Af* in 
Fig. O Related parameters used for this example are as follows: log( '^'^"^^ ) = —2, pi/p2 = 5, 
mi/m2 = 1.2, a/d = 125, Rja = 0.02, = 0.01, A = 35°, i = 89*5°, and 9 = 10°. See 

details in Section 13.21 
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The transit of a binary is different from that of a single body. As illustrated above, the binary 
may enter or exit the transit one by one, or the projected areas of the two components may 
overlap, so that some special features can be shown in the transit light curve (e.g., some step 
changes or "bulges"). Some properties of the binary planet (e.g., radii o f its two components . 



and its angular velocity and semimajor axis if its u' is sufficiently fast; ISato fc Asadal 120091 ) 
can be extracted from the features. In addition, the transit duration variation and the 
transit timing variation measured from the light curves have also been propo sed to obtain 



the exomoon mass and the semimajor axis of the moon 's orbit JKippinejboOQal jbll Below 



we 

illustrate that the orbital configuration of the binary planet can be further constrained by 
the evolution curve of the observational stellar radial velocity. 

Note that the velocity in Figure [3] (similarly in Fig. Elbelow) is expressed in a dimension- 
less quantity, where the stellar parameters involved in the normalization could be non-trivial 
to estimate in realistic systems and would be done through other independent methods and 
abundant knowledge in stellar astrophysics. 



3.2. Stellar radial velocity anomaly and orbital configuration of a binary 

planet 

As mentioned in the Introduction, a binary planet may leave signature on the observa- 
tional radial velocity of the star, as well as on its light curve. The apparent stellar radial 
velocity anomaly due to the blocking of the stellar light is given by 

./ / xl{x, z)dxdz 



J J I{x, z)dxdz 



(5) 



[e.g.. lOhta. Taruya. fc Sutdl2005uWinn et al.ll2005l ). where K = f2*sini^, is the line-of-sight 
component of the spin angular velocity and may be observationally constrained from the 
stellar spectrum. Applying equation (|5]) to the example shown in Figure [3^, we obtain the 
evolution curve of its corresponding radial velocity anomaly in Figure [Sb (see also Fig. 1 in 
Simon et allboool ). 



Below we demonstrate how the dynamical configuration of the binary planet is incor- 
porated into the evolution curve of the stellar radial velocity anomaly. For simplicity, we 
consider the complete transit phase of the binary planet. And we analyze the case in which 
the projected areas of the two planets do not overlap, and the non-overlap is likely to be 
true during most of the transit time especially ii Ri + R2 d oi \t: /2 — i'\ > (_Ri + R2)/d. 
Thus, the stellar radial velocity anomaly can be simplified as follows: 

xiAi + X2A2 



Av. = K 



A^-Ai- A2 



(6) 
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where xi = xq' + x[ and X2 = xq' + X2 are the x-coordinates of the center of each planet, 
respectively (see Fig. |2)i), 



Xo' = a(l — sin^ A sin^ i)^^"^ sin[a;(t — to)] (if e = 0) 



(7) 



is the x-coordinate of the center of mass of the binary planet, to is set so that xq' = at 
t = to, 



and 



X2 



di{l- sin^ A' sin^ i'Y^^ sm[uj'{t - to) + 4 
-^2(1 - sin^ A'sin^ i'Y^^ sm[uj'{t - to) + 



(9) 



are the x-coordinates of the two planets relative to their center of mass, di and ^2 are the 
distances of the two planets to their center of mass, and represents the orbital phase of the 
inner binary. The y-coordinates of the planets along the line of sight do not appear explicitly 
in equation ([6]), but they are involved in the expression through the angles u'{t — to) and (p. 
Applying equations ([7])-([9]) to equation IQ, we have 



Av^ = Avo' + Avb, 



where 



Avo' = K (1 - sin^ A sin^ i)^'^ sin[a;(t - to)], 

^ ^il^iiZ:ML(i_sin2A'sin2 0'/^in[a;'(t-to) + 
The Af* is composed of the two terms, Avqi and Av^. 



(10) 

(11) 
(12) 



The Avqi gives the radial velocity anomaly as if the transiting body is a single body 
with the values of its mass and projecting area being the total ones of the binary and 
contains the orbital configuration of the outer binary, i.e., the angles A and %. This 
R-M effect due to the transit of a single planet has been used to extract those angles of 
some realistic exoplanetary systems. Together with some assumption or observational 
evidence on the possible distribution of the inclination of the stellar spin i^,, the orbital 
inclination of the planet relative to the stella r spin (i.e., the an gle between n and n*) 
can be further statistically constrained (e.g., IWinn et al.l 120091 ). This inclination has 
also been measured in a number of realistic systems through the photometric anomalies 
exhibited in transit light curves, which are interp reted as passages of the planet over 



dark starspots (e.g.. ISanchis-Ojeda fc Winnll201ll ). 
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• The Aff, comes from the relative motion of the inner binary, for which the sine mode of 
equation (|T2l) represents its periodical orbital motion. The information on the orbital 
configuration of the inner binary (A', i') comes from Avf,, and we focus on the effects 
of this term in this paper. 

Within one transit duration, as u6t ~ -R*/a ^ 1, the Avq' in equation ([TT]) can be 
simplified to be linear with time as follows, 

Avo' = K + (1 - sin^ A sin^ ^flMt - U)- (13) 

A^: — /ii — /i2 

If the sine mode in Av^ can be identified in the observations, its period and amplitude can be 
used to constrain the value of u' and the geometric configuration (A', i'). If u}'5t ~ ^ 
the Avb in equation (fT2|l can be reduced to be also linear with time as follows 

= i^Aill^M^(i_sin2A'sin2^')V2[sin0 + c^'(t_to)cos0], (14) 
— Ai — A2 

and the slope of the Av^-t curve during the complete transit phase is given by 

dAv^ 

kq' + kh cos 0, (15) 



dt 
where 

= ir-^^^i-t^(l - sin^ Asin2^)V2^ (16) 
A^-Ai-A2^ ^ ^ ^ 

and 



7 f ^1^1 ~ ^2^2 A /, ■ 2 \/ • 2 •/\l/2 / 

kh = K \ — — 1 — sm A sm i) ' u . 

\A,-Ai-A2j^ 



(17) 



The koi is constant with time as A and i. The fcb depends on the orbital configuration of the 
inner binary (A', i'), which is usually constant within one transit duration, as mentioned in 
Section [2J 

Even if the eccentricity of the outer binary e is non-zero, but has a low or moderate 
value [e.g. e=0.3, so that the linear approximation of uj{t — to) in eq. [13] is still valid], the 
expressions for Avq' and ko' can be modified simply by replacing a with ga in equations 
f lTS]) and f lT^ . where the factor 

^ = (1 - e2)-^/2[l + ecos(a + ^POQ)], (18) 

the angle a is defined by tana = sinAcosi/ cosA and cosa = — cosA/(l — sin^ A sin^ i)^/^, 
and the meaning of the angle ZPOQ is indicated in Figure [TJ 

In a longer time period At{> 6t), the evolution of A' and i' is determined by the 
precession of n around n' and the value of the angle 6. In this case, we discuss the linear 
mode in the following two regimes. 
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• If QAt ^ 1 (i.e., At is much shorter than the precession timescale, but covers multiple 
transit events), ki, is still roughly constant. The distribution of the phases 0/27r follows 
the distribution of nu'/u [n: integer). In general, unless u'/u is an integer, is 
uniformly distributed between and 27r, and we can get the average of the slopes over 
the multiple transits as follows 



which can be used to constrain the sky-projected angle A. The rms of the slopes is 



The value of 6* is a function of A' and i' (see eq. [T]), and equation fl^U]) can be used as 
an observational constraint for statistical determination of the probability distribution 



• If QAt > 1, kf, changes due to the precession of n and the variation of A' and i'. Here 
we only discuss the case with < ^ < 40° or 140° ^0 ^ 180°. As seen from Figure HI 
the evolution pattern of / = (1 — sin^ A' sin^ i')^/^ is different with different 9. Note 
that the evolution pattern for vr — 6* (0 < 6* < 40°) can be obtained by reversing the 
time in the pattern for 6, as the precession direction for 6 > 90° is along —n. The slope 
of k should distribute within the envelope of /, and the magnitude of the variation of 
/ increases with increasing 6 (for 6 < 90°). Hence the distribution of the slopes can 
be used to constrain the value of 6 (or |90° — ^^|). We illustrate such an example in 
Figure |5]^c)-(e), by doing full three-body numerical simulations on dynamical evolution 
of a binary planet rotating around a star and obtaining its multiple transiting events. 

Similarly, if the stellar radial velocity anomaly due to a binary planet Avb is in a sine 
mode in equation (IT^ . observations of multiple transit events within a longer time At can 
be useful to obtain the evolution of the amplitude of the sine mode, the evolution of the 
orbital configuration of the inner binary (A', i'), and also further statistically constrain the 
angle 6 (see the example shown in Fig. |5^-b). 

In addition, in some cases, the binary planet are not in the complete transit phase, but 
only one component is in the transit. For example, this may occur in the ingress or egress 
phase; and ii d > 2i?*, it is also more likely that each component transits in front of the 
star one by one. In these cases, it is easy to generalize the analysis above and obtain the 
deviation in the R-M effect due to the transit of each component, by setting the projected 
area of the other component to zero in the formula above. 



k = ko' 



(19) 



{k - ky 



1/2 




(20) 



of ^. 



Note that in the modeling of most of the other exomoon detection methods, the orbit of 
the exomoon surrounding the planet is assumed to be co-aligned with the orbit of the planet 
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Fig. 4. — Examples of the orbital evolution of a binary planet rotating around a star, 
calculated from a full three-body numerical simulation. The related parameters used are the 
same as those in Fig. |3] except the value of 6. The variable of f{6, t) affects the amplitude or 
slope of the stellar radial velocity anomaly curve (see eqs. [T^ and [TTIl . Panel (a) is mainly 
for 6 < 40°; and panel (b) gives an example with 6 = 75°, where the Kozai mechanism 
is effective. In panel (a), different curves of / represent the evolution under different 6, 
which produces different slope distributions of the stellar radial velocity anomaly curves, 
given multiple transit observations. Panels (c) and (d) show the evolution of some other 
parameters corresponding to the cases shown in panel (a) and (b), respectively. In panel (b) 
{6 = 20°), the parameters e', 6, and Ck are almost constant with time, and the angles A' 
and i' evolve periodically due to the precession of the orbital angular momentum. In panel 
(d), the eccentricity e' and the angle 6 exchange periodically, while Ck remains constant. 
See details in Sections [2] and [221 
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Fig. 5. — The Aff, of the deviation in the R-M effect due to the transit of a binary planet 
(see eq. [12]). Panels (a)-(b) are for the sine mode of Aff,, and panels (c)-(d) for the linear 
mode. For each mode, we show the results for different 6. The figure is obtained by full 
three-body numerical simulations. The related parameters used in the simulations are the 
same as those in Fig. [31 except log( "^i+"^2 ^ _ _4 f^j, ^j-^g linear mode and the different 
values of 9. In each of panels (a)-(d), each curve represents the Av^ within one transit 
event of the binary planet, and multiple curves represent the results of its multiple transit 
events. We obtain 500 continuous transit events for each system and display only 20% of 
them (chosen randomly) in those panels for visual clarity. Note that the time span of the 
transit events is longer than the precession period 27r/f2, and a longer time span can give 
a better statistical illustration in panel (e). Each curve of the transits is different due to 
the time evolution of the angles A' and i' and/or different orbital phases (f). All the curves 
in each panel are stacked up for comparison of either their amplitudes or their slopes. In 
panels (a) and (b), the curves have been shifted along the x-axis to have Aff, = at t = 
and have their corresponding maxima/minima aligned; and in panels (c) and (d), the curves 
have been shifted along the |/-axis to have their middle points at the zero point, and the 
constant Ci on the ?/-axis labels indicates the shifts, which are different for different events. 
In each curve, only the complete transit phase is shown with the overlapping period of the 
two planets being removed for visual clarity, which can be done based on its corresponding 
transit light curve, although in reality some system information can also be extracted from 
the ingress, egress, and overlapping parts. The projection of a binary planet is more likely 
to overlap along the line of sight when 9 is small. The gap in the middle of the curves shown 
in panel (a) is just the result of removing the overlapping phase of the two planets. As seen 
from panels (a) and (b), the scatter of the velocity amplitude in the sine modes is larger for 
the larger angle 9 = 40° than that for 9 = 10°, which is consistent with the behavior of / 
shown in Figure [H Panel (e) gives the slope distributions of the Avb-t curves obtained for 
the linear mode, which also appear different for different 9. 
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surrounding the star. ISato &: Asadal (120 lOl ) consider how the special step or overlap features 
shown in transit light curves (cf., Fig. [3^) can be used to infer the orbital inclination of an 
exomoon, in which the relevant cases ar e for the condition that i' is close to 90°. The transit 
duration variations derived by iKippingl (j2009b|) involve the different inclination parameters 
of an exomoon, but which is limited to uj'St <C 1 or the linear mode and ignores the evolution 
of the dynamical systems (i.e., the precession of n'). 



3.3. Systems likely to be revealed by observations 



To have the signatures of a binary planet detectable in observations, the change of the 
stellar radial velocity anomaly due to the binary planet should be significantly large during 
each transit. As analyzed above, the deviation Aw is composed of two parts, Awo' and Afb 
(see eq. [TOj) . and both have the contributions from a second planet or exomoon. 

For the first part Afo', the contribution from each component of the binary can be 
estimated by 

KR^ \ /lOMr 



5ms ^ 



1,2, 



(21) 



5km s ^ J \ 

which is the same if the two components have identical projected areas; and the special 
features (e.g., step change EE') shown in Figure 3 have the same orders of magnitude as 
that estimate and may serve as some characteristic signals of binary planet candidates in 
the R-M effect. 

For the second part Af^, which reveals the dynamical configuration of the inner binary, 
we used some individual systems with specific dynamical parameters to indicate its effect in 
the above section. Here to see a general parameter space of binary planet systems that are 
likely to be revealed in Av^ by observations, we define the following velocity change: 



5vh 



K 



A^- Ai- Ai 



u'5t, if uj'5t < 1, 



.dJ + mi + 777-2 



1/2 



Ai + A^ 



f, 



and 



diAi - ^2^2 



A^- Ai- A2 



K 



KR.{ — 



if u'6t > 1, 



A1 + A2 

A^: 



(22) 
(23) 

(24) 
(25) 
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where 




d 



A1+A2 

A* 



1 - 



1 + 1711/7712 



1 



1 + 



( 



mi/m2 

Pl/P2 



1 



) 



2/3 



if 



A1+A2 



< 1. 



(26) 



Equation fl22l) represents the maximum change of the stellar radial velocity anomaly due to 
the inner binary in the linear mode of equation f|T^ . i.e., kbSt; and equation flMl) represents 
the amplitude of the sine mode shown in equation f[T^ . Note that the 6vb in both equations 
fl22]) and is a defined variable. Although the expressions are obtained through the limits 
at ijj'6t <^ 1 and u'6t ^ 1, their values at (x)'St should be in the transition of the two limits 
and the equations above work for the order-of-magnitude estimates and the purpose of the 
paper. As seen from the equations above, not only the planet areas/sizes are involved in 
the R-M effect as for a single planet, but their masses or mass densities are also involved 
in the effect for a binary planet due to the relative motion of the two components of the 
binary, as the hidden areas of the stellar surface are affected by the relative positions of 
the binary components and the relative distance of each component from the center of mass 
of the binary (cf., eqs. |8] and |9]) is determined by its two component masses. Given the 
mass ratios of the two components, the radius ratios can also be expressed through their 
mass density ratios. Thus, the amplitudes of the R-M effects indicated in Figures EHS] are 
expressed through the extra dimensions in multiple panels. The involvement of the extra 
(mass density) dimension in the study is useful especially considering that recent Kepler 
observations have revealed the mass density of planets do span a large range. In addition, 
the R-M effect for the transit of a single planet is related with the orbital semimajor axis, 
but the semimajor axis a of the outer binary is involved in the effect for a binary planet 
as shown in equation f l23|) because the transit time 5t is affected by a and a longer transit 
time leads to a larger change of the relative position of the binary components and further 
a larger deviation in the stellar radial velocity. 

As seen from equation fl26l) . the value of 6vb is significant only if the two components 
of the binary planets are different, especially if the heavy one has a high mass density and 
the light one has a low mass density. If the two components have the same mass and mass 
density, they have the same projecting area and their motion is symmetric, and thus we have 
6vb = with fsv^ = 0, although in this case the contribution from Avq' (eq. E]) can be large 
due to the combination of the projected areas of the two single components. In addition, the 
value of 6vb is significant only if the area of the binary planet Ai + A2 is significantly large, 
as 5vb is proportional to ^^^^^ . 
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By using equations ( l23l) and fl25|l . we show how Svb depends on other various properties 
of the system in Figures [6H71 We set the parameters KR^ = 5km s~^ and ^^^f^ = 0.01 
in both of the figures. We set -R*/a = 0.005 and 0.05 in Figures [H] and [TJ respectively. If 
R^: = Rq, R^,/a = 0.005 corresponds to the distance of the earth from the sun (a = lAU), 
and i?*/a = 0.05 corresponds to some typical semimajor axis of hot Jupiters discovered in 
the vicinity of a star (a = O.lAU). The contours of 6vb as a function of (mi + 1712) /m^ and 
a/d are shown by black solid curves. In each figure we display how the contours change 
for binary planets with different mass ratios {mi/ 1112 = 1.2, 10, 100, 1000) and different mass 
density ratios (pi/p2 = 5,1,0.2). As seen from the figures, in the region below the black 
dotted curve (i.e., u'St < 1), 6vb increases with increasing (mi + m2)/m* and a/d, as u'6t 
does in the same tendency (see eg. [23|l : and in the region above the curve, 6vb increases with 
decreasing a/d (given R^/a), as the difference of the line-of-sight rotation velocity of a star 
covered by each component is likely to be relatively large for a wide binary planet (with large 
d). By comparing the contours obtained for different mass density ratios, the figures also 
illustrate that the values of Svb can be significant (e.g., up to ms~^ or several ten ms~^) only 
if the two components of the binary planets are different, especially for high pi/ p2 ratios (see 
panels la-lc), as analyzed above. Given pi/p2, the 6vb at bottom panels (e.g.. Id, 2d, 3d) 
have relatively low values, which indicates the difficulty to detect the effect of a too small 
exomoon. Figure [7]has a higher R^/a than Figure [61 On the one hand, the curve of uj'6t = 1 
in Figure [7] (black dotted curve) shifts downwards; and thus, although the contours of 6vb 
below the black dotted curve are the same in both of the figures, the region above the curve 
has relatively low 6vb in Figure [71 On the other hand, a higher R^,/a would imply a shorter 
period of the binary rotating around the star and thus a higher probability to do multiple 
transit observations to get a better statistics for the system. In addition, for different values 
of KR^, and the contour values of 5vb in the figures should be adjusted simply by 

multiplying them by a factor of { ^^Jlu ) ■ i^^^)- 

Based on the results above, we discuss the magnitude of the 6vb in the following examples 
of binary planet or exomoon systems. We assume that the central star is a solar-like star 
(with solar mass and radius). 



An earth-moon system (with mi/m2 — 81, Pi/p2 — 1-6, R^/a ~ 0.005, a/d ~ 400) 
may cause a 5vb only up to 0.1cm s~^, e.g., due to the small sizes of the earth and the 



10 



-4^ 



moon (with ^'^^^'-^ 
to the estimate by equation ([21 
can be up to 1cm s~^ (consistent with the value shown in [Simon et al.ll2010l . where the 
effect of Avb is not discussed). 



which is too small to be detected. Note that according 
, the contribution of th e moon to the dev iation Av'q 



The Ganymede is the biggest moon in the Solar system. A Jupiter- Ganymede system 
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(with mi/m2 — 1.3 x 10^, pi/p2 — 0.7) may lead to a 5vh only up to several cms ^ 
(for Rja - 0.001, a/d ~ 800). 

• For a Jupiter-rocky moon system (e.g., with pxj pi = 0.2): a moon or satellite with 
mass 1712 < O.lmi has a low 6vb less than ms~^ (cf., panels 3b-3d in Figs. |6H7]). A 
Jupiter-earth system (with m2 ^ 3 x lO^^rrii) has a low 6vi, only in the range of 1- 
10cm s~^ (see panels 3c-3d). However, if the rocky moon/satellite has a mass close to 
the Jupiter (see panel 3a), the 6vb can be up to several to several ten ms~^, which 
may be detectable by the current techniques. Note that this is an exotic all 
the rock bodies revealed by the Kepler do not exceed several ten earth mass so fafl 

• A binary Jupiter-like planet system (e.g., with component masses mi = lOMj and 
7712 = 1 or O.lMj, and with Jupiter-like mass densities pi = P2) may have a 6vb 
ranging from ms~^ to several ten ms~^, although hot binary Jupiter systems (with 
relatively high Rt^/a) have relatively low 5vb- The 6vb can be even larger if pi > p2- 
For such a system, if any, its signatures on the R-M effect may be detected by future 
observations. 

• The CoRoT-9 system has a solar-like central star and its orbiting exoplanet CoRoT-9b 
has the mass and radius close to the Jupiter's. The CoRoT-9b is one of the longest 
period transiting Jupiter (~ 95 days) that has so far been confirmed and has a semi- 
major axis a ~ Its i?*/a ~ 0.01 is between the cases shown in Figs. Eland 
[71 As inferred from the figures, if the CoRoT-9b has a satellite m2 > 10~^mi, its 
5vh can range from several cms~^ to several ms~^, depending on the detailed satellite 
properties. 



By using equations ( !T3|) . (fT6|) . fl22|) . and ( 12^ . the ratio of the two parts in the deviation 
Af * (see eq. [T0|) is about 

^Vb Q ^ / 6vb \ /5kms^\ / O.OIA, \ ^^7) 



ko'dt V5ms-V V KR^ J V^i + -4 

As mentioned above, both of the two parts have the contribution from a second planet or 
exomoon, and the ratio of the two contributions can be estimated by ( sm" ) ( ^ kt^^ ^ ) ( ^° a^* ) ■ 



•^http://kepler. nasa.gov/ 
"'http://exoplanet.eu/ 
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Fig. 6. — Equi-5vb curves for different '"^^'"^ and a/d. The contours of 6vb are drawn in 
black solid curves and calculated from equations ( 12^ and (^B^. The values of Svb labeled for 
each curve are in units of ms~^, and they are chosen in an interval of A log(5f6/ms~^) = 0.3. 
The cyan dotted curve indicates the Hill radius (eq. [2]), and a binary planet with parameters 
located below the curve cannot survive due to the tidal breakup by the star. The black 
dotted curve indicates u'6t = 1, which separates the two regimes expressed by equations 
fl25]l and (1251) . The region below the black dotted curve has u'6t < 1. The parameters are 
set as follows: KR^: = 5km s~^, "^^^^^ = 0.01, and = 0.005. The red line is a reference 
line for d = 2R^. The Sv^ labeled are estimated for a binary in the complete transit phase. 
Generally a binary located to the left of the red line is more likely not to have the complete 
transit phase and the transit of each component in front of the star occurs one by one; in 
this case, the 5vb can also be easily estimated, and one component may have a larger value 
than the labeled one and the other component has a smaller value. As labeled in the figure, 
different panels give different mass and mass density ratios of a binary planet, i.e., pi/p2 = 5, 
1, 0.2 and mi/ma = 1.2, 10, 100, 1000. 
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Fig. 7. — Same as in Figure El except that the parameter ^ = 0.05 is set in this figure. The 
cyan dashed fine shown in some panels represents the Roche hmit, calculated from equation 
([3]); and the parameter space of a binary planet should lie to the left of the line. The blue 
dotted line is a reference line for d = Ri + R2, and a binary planet should also lie to the left 
of the line. As in Figure |6l the grey shaded areas indicate the parameter space that a binary 
planet cannot survive dynamically. Both the cyan dashed line and the blue dotted line are 
not shown in Figure [HI as they he beyond the upper bound of the x-axis of each panel. 
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Note that a binary planet has a different gravitational effect on the stellar motion from 
a single planet. The quadrupole moment of the gravitational force from the binary planet is 

and its effect on the dynamical motion over each transit duration can be estimated by 

^ . loms-' (^^) f 12^V f W , (29) 

VlOOms-i/ \ a J \ mi J \ a J ' ^ ^ 

which is generally negligible. The v^, is the velocity of the star relative to the center of mass 
of the system. It would be interesting to investigate the long-term dynamical effect of the 
binary planet on the stellar radial velocity, but which is beyond the scope of this paper. 



3.4. Application to hierarchical triple star systems 

Figures [MZI can be applied to a hierarchical triple star system (e.g., for the parameter 
space '"^^t"™^ ^, 10~^i in which a dark binary star is transiting in front of a tertiary star 



(e.g., 



Carter et al.ll201l[ ). Here by "dark" we mean that the light emission from the binary 
is ignored as planets for simplicity. It is easy to generalize the analysis above to include 
the light emission from the binary. If the transiting binary is a compact object (e.g., white 
dwarf or neutron star) plus a planet /brown dwarf, the factor fsv^ is up to 1, and the Svb can 
be much larger than that of binary planet system. Figure |8] illustrates such a case. As seen 
from the figure, the Svb is high, up to lO^ms"^. Thus, such hierarchical triple star systems 
may be revealed through the R-M effects by future observations. 

Figure [9] illustrates two examples of the transit light curves and the radial velocity 
anomaly curves for the triple star systems with dynamical parameters located in the param- 
eter space shown in Figure [H As seen from the figure, the radial velocity anomaly caused by 
the rotational motion of the planet/brown dwarf around the compact object is much more 
significantly displayed either through the bulge/hill/trough features during the ingress/egress 
phase (for panel a) or in the sine-like curve during the complete transit phase. The values 
of Avb can be comparable to Avq' ■ It is plausible to expect that the features of these curves 
are very useful to infer the dynamical configurations of the systems, as illustrated for binary 
planets above. A more detailed discussion on extracting the configurations is beyond the 
scope of the paper. 
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Fig. 8. — Equi-Svb curves for hierarchical triple star systems. The curves have the same 
meanings as those in Figures [6] and [71 The values of 5f f, labeled for each curve are in an 
interval of Alog((5vb/ms~^) = 0.15. The binary star in the hierarchical triple star system is 
composed of a compact object (e.g., white dwarf or neutron star) plus a planet/brown dwarf, 
so a very high density ratio pi/ p2 is used in the figure. The parameters KR^ and ^^^^'^ 
are the same as those in Figure [61 The i?*/a is 0.005 for panel (a) and 1/6 for panel (b), 
respectively. As mxjm^ decreases, the cyan dashed line shifts rightwards, and the contours 
of bvh are affected little for mi/m2 > 10. See also Section [3^ 
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Fig. 9. — Examples of the simulated transit light curves and the stellar radial velocity 
anomaly curves for a dark binary star transiting in front of a tertiary star. The methods 
to obtain the simulated curves are the same as those for Figure [31 The system parameters 
used in panels (a) and (b) are chosen from the parameter space shown in Figure [Ht^a) and 
(b), respectively. Panel (a) has log(^^^^) = 0, a/d = 150, i = 90°, and 9 = 10°; and 
panel (b) has log(^^^^^) = 0, a/d = 4, i = 89.5°, and = 10°. Both of the examples 
indicate a significant amplitude of changes in the radial velocity anomaly caused by the 
relative motion between the compact object and the planet/brown dwarf. In panel (a), a 
number of the bulge/hill/trough features are indicated in the curves, which is because the 
orbital period of the dark binary is relatively short so that the planet/brown dwarf rotating 
around the compact star can move completely into, and then move out of, and then re-move 
into the projected stellar surface before the complete transit phase. The similar effect of the 
rotational motion also occurs during the egress phase. In panel (b), the sine mode of Avb 
represented by equation ([T2|) is significantly shown in Av^. See also Section [331 
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4. Discussion 

The properties of binary planets (e.g., mass, size, and semimajor axis) can be constrained 
through their signatures on the transit light curves and stellar radial velocity curves. Our 
analysis of the R-M effect for a binary planet during its complete transit phase show that 
effect is composed of two parts. The first part is the conventional one similar to the R-M 
effect from the transit of a single planet with the combined masses and projected areas of 
the binary components (eq. [11]); and the second part is caused by the orbital rotation of 
the binary components, which may add a sine- or linear-mode deviation to the stellar radial 
velocity curve (eq. [12]) . 

In this paper we focus on the discussion on the second part, and we find both the 
evolution of the orbital rotating phases and the precession of the binary orbital plane may 
lead to different amplitudes of the deviations in different transit events of the same system. 
The resulted distribution and dispersion of the deviations in multiple transit events can be 
used to extract the orbital configuration of the binary planet or even the inclination of its 
orbital plane relative to the plane of its center of mass rotating around the star (e.g.. Fig. [5]). 

The second part of the R-M effect is more likely to be revealed if the binary components 
have different masses and mass densities, especially if the heavy one has a high mass density 
and the light one has a low density. For example, our calculations show that the signature 
can be up to several or several ten ms~^ with the mass ratio mi/m2 up to 10^, if the mass 
density ratio pi/p2 = 5. A small and rocky exomoon with m2 < O.lmi would cause a low 6vb 
less than ms^^. A strong signature may be caused if at least one of the components of the 
binary planet is a giant planet. Note that stellar noises produced by oscillations, granulation 
phenomena, and activities coul d contribute to the cha nge of stellar radial velocities with 



an amplitude up to ms ^ (e.g., jPumusque et al.ll201ll ). but they have their own variation 



periods and patterns to be distinguished from the effect of binary planets, and some statistical 
methods can be developed to extract smaller signals from the noises. 

A long observation time would cover multiple transits of a binary planet, and a better 
statistics on the distribution of the deviation in the R-M effect could be potentially obtained. 
To trace the evolution of the geometric configuration of the system, we need at least one 
period of the orbital angular momentum of the inner binary precessing around that of the 
outer binary. The precession period is roughly about u'/u times the orbital period of the 
outer binary. For the parameter space shown in Figures [HI and [7] u'/u is generally less than 
200 (cf., the oj'St = 1 curve, which can also be taken as a reference line for u'/u ~ {R^/a)^^). 

The misalignment between the plane of the binary planet and its rotating plane around 
the star is one of fundamental parameters of the dynamical system. A large misalignment is 
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likely to cause a large dispersion or different distribution of the deviations in the R-M effect. 
Recent measurements have discovered that the orbit of a planet may be highly inclined to the 
stellar spin. Different mechanisms have been proposed for the formation of those misaligned 
orbits , e.g., Kozai capture, planet-plan e t scattering, resonance capture by planet migration 



Murrav-C 



Yu fc Tremaine 



av &: Schlichtingl 1201 It iNagasawa et al.l l2008t iFabrycky &: Tremaine 



2007 



200ll ). Similarly, the different configurations/inclinations of binary planets 



(see 6 defined in Tabled]), if detected in future, should be also useful in constraining forma- 
tion mechanisms of binary planets or exomoons, and shed new light on our understanding 
the diversity of planetary systems. For example, different moon formation mechanisms 
may lead to their different kinematic distributions. Moons formed from the disk material 
surrounding a planet are predicted to have prograde orbits, while those formed from gravita- 
tional capture /impacts /exchan ge interactions can have either prograde or retrograde orbits 
(jjewitt &: Haghighipourl I2OO7I : and references therein). It is also likely that the orbit of a 
moon could be affected by the later evolution of the system (e.g., by the later inner/outer 
migration of the outer/inner planets). Regarding a binary planet syst em with comparable 



comp onent masses, the study of their formation theory is starting (e.g., iPodsiadlowski et al. 



2OIOI ). and one of the most exciting steps would be to discover a realistic system in observa- 
tions in the near future. 

To discover binary planets and exomoons becomes promising and practical with future 
developments in instruments, which also lay the foundation for finding the signatures dis- 
cussed in this paper. For example, planned ground-based surveys such as the Large Synoptic 
Survey Telescope may detect thousands to tens of thousands of planetary transit candidates; 
and the space missions, PLAnetary Transits and Oscillations of stars and Transiting Exo- 
planet Survey Satellite, aim to find transiting planets around relatively bright stars, making it 
easier to confirm discoveries using follow-up radial velocity measurements. Hopefully follow- 
up observations could provide some binary planet or exomoon candidates. The astro-comb 
technique is air ning to achieve a precision as h i gh as 1cm s~^ in astronomical radial velocity 
measurements (ILi et al.ll2008l ). iKipping et al.l (|2012| ) also proposed a systematic search for 
exomoons around transiting exoplanet candidates observed by the Kepler mission. 

We studied the basic signatures of binary planets that are likely to be revealed in the 
R-M effect. To understand the roles of the crucial parameters played in the signatures, 
we have made some approximations in our study, which could be improved or adapted to 
reahstic systems in future work, for example, the study could be extended to a general case 
in which the center of mass of the binary planet is on an eccentric orbit. The study would 
become complicated if an exoplanet has multiple moons. In this case, the small moons would 
contribute little to the deviation in the R-M effect, and the one with a relatively large radius 
and located at a relatively far distance from the primary planet would imprint the most 
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significant effe ct. The limb- darke ning effect can reduce the amphtude of the R-M effect by 



20-40 percent (ISimon et aLll2010l ). and also affect the shape of the radial velocity anomaly 
(for both the linear and the sine cases studied in this paper). This effect could be corrected 
with the aid of the observational transit light curves in reality. A statistical method to map 
the reconstruction of relevant parameters of a binary planet system is beyond the scope of 
this paper, but would need to be explored in details in future. 

After extending the results to a hierarchical triple star system containing a dark binary 
and a tertiary star, the deviation in the R-M effect would be large enough to be detected 
especially if the dark binary is composed of a compact object and a brown dwarf/planet, 
which may put further constraints on the geometrical configuration of triple star systems 
and provide insights on their formation and evolution. 

We thank the referee for many helpful comments. This research was supported in part 
by the National Natural Science Foundation of China under No. 10973001. 
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